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Solvent structure and dynamics of acetonitrile at its liquid-vapor (LV) interface and at the acetonitrile-silica 
(LS) interface are studied by means of molecular dynamics simulations. We set up the interfacial system and 
treat the long-ranged electrostatics carefully to obtain both stable LV and LS interfaces within the same system. 
Single molecule (singlet) and correlated density orientational profiles and singlet and collective reorientational 
dynamics are reported in detail for both interfaces. At the LS interface acetonitrile forms layers. The clos- 
est sublayer is dominated by nitrogen atoms bonding to the hydrogen sites of the silica surface. The singlet 
molecular reorientation is strongly hindered when close to the silica surface, but at the LV interface it relaxes 
much faster than in the bulk. We find that antiparallel correlations between acetonitrile molecules at the LV 
interface are even stronger than in the bulk liquid phase. This strong antiparallel correlation disappears at the LS 
interface. The collective reorientational relaxation of the first layer acetonitrile is much faster than the singlet 
reorientational relaxation but it is still slower than in the bulk. These results are interpreted with reference to a 
variety of experiments recently carried out. 

In addition, we found that defining interface properties based on the distribution of positions of different 
choices of atoms or sites within the molecule leads to apparently different orientational profiles, especially at 
the LV interface. We provide a general formulation showing that this ambiguity arises when the size of the 
molecule is comparable to the interfacial width and is particularly significant when there is a large difference 
in density at the upper and lower boundaries of the interface. We finally analyze the effect of the long ranged 
part of the electrostatics to show the necessity of properly treating long-ranged electrostatics for simulations of 
interfacial systems. 



I. INTRODUCTION 

Interfacial properties of liquids has been a challenging and 
important subject to study because both structure and dynam- 
ics of molecules at interfaces are usually quite different from 
those in bulk liquids. In the past two decades, several exper- 
imental groups have focused on the acetonitrile system con- 
fined in silica nanopores or at its air/liquid interface using a 
variety of experimental tools iflT- fToh . Early NMR experiments 
by Zhang and Jonas suggested a two-state model, which says 
that the liquids confined in pores with diameter 24 A or 44 A 
consist of both an adsorbed component and a bulk-like com- 
ponent yfl. Additional support for this picture was provided 
by Fourkas and coworkers, who employed optical Kerr effect 
(OKE) spectroscopy to study the dynamics of liquids confined 
in silica nanopores. Their experiments suggested that hydro- 
gen bonding is responsible for the extreme inhibition of dy- 
namics at the pore surface JH-Q]. Some shortcomings of the 
two-state model in describing detailed features of the dynam- 
ics have been noted and addressed iffl fllll . 

Other experimental work using isotherm adsorption tech- 
niques has been explained by dividing liquid acetonitrile 
molecules confined in pores into monolayer acetonitrile and 
capillary condensed acetonitrile J3,Ht]. Kittaka and coworkers 
have recently carried out both spectroscopic and temperature- 
gradient experiments to determine low temperature properties 
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of acetonitrile confined in mesoporous material MCM-41 @]. 
Their measurements have shown that the monolayer acetoni- 
trile adsorbed onto the inner surface of the cylindrical pores 
interacts strongly with surface hydroxyls, in agreement with 
the earlier finding from the OKE experiment. Their work and 
the earlier Raman experimental work by Tanaka and cowork- 
ers J3l both pointed out that the rotational relaxation time 
scale of the remaining capillary condensed molecules inside 
the pores is similar to that of the bulk liquid phase. They have 
also found that the melting point in the confined silica pores 
from phase transition of a-phase crystal to liquid acetonitrile 
is much lower than the value measured in the bulk. 

Surface sensitive experiments have also examined the 
liquid-vapor (LV) interfacial properties of acetonitrile and the 
binary mixture of water and acetonitrile 12j,|8[]. These experi- 
ments aim to gain information on molecular ordering at the in- 
terface where the molecular environment is highly inhomoge- 
neous. Using the sum frequency generation technique, Eisen- 
thal and coworkers showed that the vibrational frequency of 
acetonitrile changes in mixed solvents as its concentration 
at the LV interface is varied JUl. Their work indicates that 
acetonitrile molecules prefer to lie parallel to the interface at 
higher concentrations. The gradual change of preferred orien- 
tation was later confirmed by Somorjai and coworkers [H but 
there is some disagreement between the results of these two 
studies J2, [Ml concerning the magnitude of the tilting angle 
and its dependence on the bulk concentration. 

A few theoretical and computational studies have been car- 
ried out to provide a molecular description of the interfacial 
structure and dynamics for both the LV and the liquid-silica 
(LS) interfacial acetonitrile systems llYlUTstl . Very recently, 
Thompson and coworkers carried out molecular dynamics 
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(MD) simulations of acetonitrile confined in silica pores with 
varying hydrobocity and analyzed in detail the structure and 
dynamics of the confined acetonitrile. Their work has shown 
that confinement and modifications in the hydrogen bond net- 
work can lead to line shape broadening and slowdown of the 
rotational dynamics JHO. A detailed MD study of the LV 
interface using both polarizable and nonpolarizable models of 
acetonitrile was also carried out recently by Paul and Chan- 
dra H3I1 . They concluded that acetonitrile molecules at the 
interface prefer to orient with their dipoles parallel to the sur- 
face. Their work also indicated that the relaxation time scale 
of molecular rotation at the LV interface is much shorter than 
that in bulk solution. 

Although both simulations and experiments on acetonitrile 
in the silica pores agree that strong surface-molecule inter- 
actions tends to slow down the dynamics of acetonitrile, it 
is not known how much of this effect is due to the cylindri- 
cal confinement and how properties would change as a func- 
tion of distance from planar interfaces. While single molecule 
(singlet) density, singlet orientational profiles and singlet ro- 
tational dynamics are often analyzed in simulation studies of 
interfacial systems, less attention has ben paid to the structural 
correlations between pairs of molecules and to collective dy- 
namics. Thus several questions remain open. How does the 
planar silica surface affect liquid correlation as a function of 
distance from the surface and how do the correlations near the 
silica surface, in the LV interface, and in the bulk differ from 
one another? Are the collective reorientational time scales af- 
fected by the inhomogeneous environment? Do they differ 
from the corresponding singlet reorientations? 

In order to answer the above questions and to provide 
a molecular description of structure, correlation, and sin- 
glet and collective dynamics, we investigate a planar sil- 
ica/liquid/vapor acetonitrile system using classical MD simu- 
lations. The potential of the planar silica surface is taken from 
the work by Lee and Rossky 1 1611 . It has also been used in 
recent simulations of water confined by silica planar surfaces 
by Giovambattista, Rossky and Debenedetti |[l7H20ll . 

To better compare our simulation with ongoing experiments 
of acetonitrile liquid on a planar silica surface 12111 . we at- 
tempt to make the simulation setup as close as possible to the 
experimental conditions. Our simulation differs from many 
other MD simulations of interfacial systems in the follow- 
ing ways: (i) Continuum implicit external potentials are used 
account for the long-ranged dispersion interaction from the 
semi-infinite bulk silica region, (ii) Long-ranged electrostat- 
ics in this slab geometry are treated with the computationally 
efficient corrected Ewald3D method originally developed by 
Yeh and Berkowitz ll22ll to correct for spurious periodic im- 
ages normal to the slab walls. The accuracy of this method is 
compared with the rigorous Ewald2D method l23ll for the di- 
mensions of our system, (iii) Our system setup allows interfa- 
cial, bulk and vapor acetonitrile to coexist in one simulation. 
Acetonitrile molecules in the system on average feel nearly 
zero pressure from the boundaries and interfacial molecules 
can readily exchange with the bulk region. 

The rest of the paper is organized as follows. Details of the 
surface construction, a validity check of the potentials, and 



simulation parameters are provided in section|II] In sectionlllll 
we investigate density distributions, orientational profiles, ra- 
dial angular correlations, and singlet and collective reorien- 
tations in detail to give explanations for a variety of experi- 
mental results. The electrostatic potential and its long-ranged 
effects are also analyzed in this section to show their relative 
importance for both interfaces. We finally draw our conclu- 
sions in section HVl 



II. SIMULATION DETAILS 

A. Surface potential and long-ranged electrostatics 

We constructed a four-layer silica surface from an idealized 
/3-Cristobalite (C9) crystal ||24ll . Silicon atoms always occupy 
the center of a perfect tetrahedron formed by the closest oxy- 
gen atoms above and below. In a real experiment, a thin sil- 
ica surface is effectively infinite on a microscopic scale. The 
positions of the atoms in the semi-infinite silica surface re- 
gion can be obtained through periodic boundary conditions 
in the xy directions and through extension of the tetrahedron 
structure in the negative z direction. The length of the x and 
y directions of the planar surface used in the simulation is 
L x = 45.605 A and L y = 43.883 A respectively and the num- 
ber of atoms in each layer of the surface is 90, 90, 270, and 90. 
The top oxygen atoms form a close-packed FCC (111) struc- 
ture with the closest O-O distance of 2a = 4-^/2/3 6, where 
b = 1.5515A is the Si-O bond length. 

We first focus on the neutral silica surface. The Lennard- 
Jones (LJ) parameters of Si and O atoms in our silica model 
were taken from the work of Rossky and coworkers lfl6[[l7ll . 
(The hydroxylated silica surface also has partial charges on 
the three atoms of the H-O-Si group and will be discussed 
later.) We used the six-site acetonitrile model of Nikitin 
and Lyubartsev where intermolecular interactions from each 
atomic site are represented by a combination Lennard- Jones 
and electrostatic potentials [25]. The usual combination rules, 
o"i2 = (ci + °2)/2 and €12 = ^(1^2 are used in the calcu- 
lation of mixed LJ interactions. By treating the semi-infinite 
silica region from the ?i-th layer to negative infinity as a con- 
tinuum, we obtain an integrated LJ 9-3 potential: 

a 1 \ — P i(J i 2?re » 2 ( ai V ( a i X 

t J n 2 VV^ a 3 [l5\*-*J \Z-ZiJ 

(1) 

where pi is the number density of atoms per unit area for the 
i-th layer, cr, and are the usual LJ 12-6 parameters for the 
type of atom in the i-th layer, and Zi is the location of the i-th 
layer in the z direction. 

The total VdW potential at position r from the whole silica 
surface including both n layers of explicit atoms and the semi- 
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FIG. 1 : The potential of a test LJ atom with parameters of the N in 
the acetonitrile model placed at four different typical sites (r a , rt, r c 
and rd) on the silica surface. (A) r a is on the top of an oxygen atom. 
(B) r-f, is on the top of the middle of two closest oxygen atoms. (C) 
r c is in the middle of the corresponding r a and rb. (D) rd is on the 
top of the center of an equilateral triangle formed by three closest 
oxygen atoms. Label for the black line is 4>c(z) = 4> c (n — 1,2). 
See Eq. © 



infinite continuum region is: 
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where N t is the number of explicit atoms of the i-th layer 
(including periodic images up to a specified cutoff distance) 
and Tj is the position of the j-th atom of the z-th layer. A 
plot of 4> n (r) with n = 0, 2, 4, and 8 respectively for a test LJ 
atom with the LJ parameters of N in the acetonitrile model as 
a function of its distance from the silica surface is shown in 
Fig. ID 

Figure Q] shows that the potential calculated using four ex- 
plicit layers overlaps that of eight explicit layers for all posi- 
tions tested. Clearly, using four explicit layers with a contin- 
uum potential is accurate enough to reproduce the VdW po- 
tential of the semi-infinite region. Note that noncharged atoms 
prefer to dock into the position r f / rather than the other three 
because of its lower energy minimum. 

The hydrogens in the hydroxylated silica surface used in 
this work are all attached to the oxygen atoms of the top sur- 
face layer (see Fig. [2] (B)). The entire silica surface (all Si 
atoms and O atoms) is fixed, with the top oxygen layer at 
the z = plane, except that hydrogens can freely rotate. 
The length of O-H bond is constrained to be 1.0 A and the 
Si-O-H harmonic angular bending potential has an equilib- 
rium angle of 109.27 degrees and a large bending strength of 
200.0kcal/mol. 




*"> W\ *""k 

l Y^4 A i i 4'V4* 

VVYYYYYV^ 

*> W2> *> Z> 3 F» 



(B) 



.•.^ O 



(D) 



FIG. 2: Side view of the simulation setup (A), with a top view of 
the hydroxylated silica surface (B). Top views of the CN (C) and 
CH3 groups (D) of the first sublayer of acetonitrile molecules dock- 
ing into the hydroxylated silica surface in a typical configuration as 
visualized by VMD software j26h are shown. 



The hydroxylated silica surface has partial charges on the 
three atoms of the H-O-Si group, which are also taken from 
the work of Rossky and coworkers lfl6l[T7ll . The long-ranged 
electrostatic interaction has to be treated carefully. We have 
carried out a test on a pair of ions in a slab-like box with the 
same length scales used in our simulation of the acetonitrile 
system (see below) using the exact Ewald2D method |p23ll . 
the Ewald3D method, and the corrected Ewald3D method in- 
troduced by Yeh and Berkowitz I22I1 . Fig. [5] shows a com- 
parison of these three methods for the test case. Clearly the 
standard Ewald3D method is very inaccurate but the corrected 
Ewald3D is in close agreement with the Ewald2D method for 
separations less than about 100 A. In what follows we use the 
corrected Ewald3D method, which is less computationally ex- 
pensive but quite accurate in the region from z = to 75A. 



B. Acetonitrile model and simulation parameters 

We performed molecular dynamics simulations of acetoni- 
trile confined in the z-direction by two parallel walls using the 
software DL_POLY version 2.18 l27ll with modifications to 
include external potentials and the corrected Ewald3D method 
to treat the electrostatics. The acetonitrile system studied here 
consisted of 864 molecules. All acetonitrile molecules were 
confined by a hydrophobic wall located at 75A and the sil- 
ica wall with its top oxygen layer located at OA. The potential 
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FIG. 3: The electrostatic interaction between two unit point charges 
computed by three different methods. The simulation setup has a 
length scale of l z = 150 A and l x = ly = 43 A.This test case is 
similar to Figure 4 in reference l22ll but for the specific length scales 
used in our simulation. 



of the hydrophobic wall had the same form as in Eq. (fTJ but 
uses only the repulsive part of the potential UM . Periodic 
boundary conditions were employed with L x = 45.605A, 
L y = 43.883A, and L z = 150A. The empty region between 
75 A and 150 A is required for the application of the corrected 
EwakBD method (see Fig. [3}. 

The simulation setup is shown in Fig.[2B). The acetonitrile 
potential energy parameters were those previously published 
by Nikitin and Lyubartsev l25ll . The cutoff distance for the 
VdW and the short-ranged part of the corrected Ewald simu- 
lation was 15 A. The screening parameter for the Ewald sum- 
mation was 0.26A . The number of fc-space vectors for the 
Ewald summation was 15, 15, and 45 for the x, y, and z direc- 
tions respectively. These parameters were optimized choices 
yielding the best efficiency for our system setup. 

A bulk simulation was initially equilibrated at T = 298K 
and P = 1 atm before putting the acetonitrile on the silica 
surface. The interfacial system was equilibrated for nearly 
a nanosecond at constant number of molecules, constant vol- 
ume and constant temperature T=298K (NVT ensemble) us- 
ing the Berendsen method rf29ll until no further energy drift 
was observed. To reduce the computational cost, we first em- 
ployed a more efficient but less accurate method ll30[[3lll with 
truncated Coulomb interactions for several hundred picosec- 
onds and then applied the accurate corrected Ewald3D method 
in the later stage of the equilibration. The latter was also used 
in all production runs. One NVE production simulation was 
run for a duration of 300ps using a time step of 1 fs. This was 
followed by another NVE simulation of 50ps using a time step 
of 0.5 fs to permit the accurate determination of time correla- 
tion functions at very short times. Positions of atoms were 
recorded at every 5 fs and every 1 fs respectively for further 
analysis. Total energy is well conserved and no temperature 
drift was observed in the NVE production simulation. 

A separate simulation of the bulk acetonitrile system was 
carried out for comparison. The production simulation of the 



bulk system had 864 molecules in a cubic box of size 42.37A. 
The number density was p B = 0.774g/cm 3 , which is in good 
agreement with the experimental value of p = 0.777g/cm 3 
at the temperature T=298K l32Tl . A top view of the surface 
and the first sublayer acetonitrile molecules on the surface is 
shown in Fig. 0(B), (C) and (D). 



IH. RESULTS AND DISCUSSIONS 
A. Singlet densities and orientational profiles 

Stable LV and LS interfaces are formed in the simulation 
of the acetonitrile system and are most simply described by 
the singlet density profile. However with larger molecules 
like acetonitrile there are several reasonable choices for which 
atomic type or molecular position should be used to character- 
ize the location of the molecule in the density profile. 

In Fig. |4]we show density profiles defined from positions 
of the nitrogen (N), methyl carbon (CT), hydrogen (H), cen- 
ter of mass (cm), and minimum z position of any atom in 
the molecule. We also defined an "identical atom" profile in 
which all atoms in the molecule are treated as if they were of 
the same type and we report the spatial distribution of this sin- 
gle atomic species. The average of the number density p(z) 
over the region 25A < z < 30A (0.98,0 s ) is nearly the same 
as the bulk density using all definitions and we consider this 
the bulk region of our interfacial system. 

Following previous work lfl3tl . the LV interfacial region is 
defined as the region over which the number density decreases 
from 90% to 10% of the bulk density. As shown in Fig. 0(A) 
and (B), density profiles from all definitions except the min- 
imum position yield almost identical LV interfacial regions 
(35.24A < z < 42.08A), while the minimum position shifts 
the LV interfacial region to 34.20A < z < 40.94A. The 
thickness of the LV interface is about 6.7A, which is slightly 
larger than the value 4.7A reported by Paul and Chandra II 13fl 
using a three-site model of acetonitrile. 

As would be expected, the LS density profiles are much 
more structured than those of the LV interface and have dif- 
ferent forms depending on the choice of defining molecular 
position. Clearly, acetonitrile forms layers on the silica sur- 
face from z = to about z ~ 20A and the various profiles 
provide information about different features of this layering. 
It is interesting that the indications of the surface-induced lay- 
ering persist past 20 A from the substrate. This contrasts with 
water near silica where perturbations die off by about 10 A 
lfl6ll . probably indicating the greater orientational flexibility 
of the tetrahedral hydrogen-bond network. 

Most density profiles show shoulder peaks except the den- 
sity profile based on the minimum position. The first layer of 
acetonitrile adsorbed on silica is defined as the region from 
z = to the first major depletion of density {z = 3A for 
the minimum position and z ~ 4.5A for the identical atom 
profile and the center of mass profile). The first minimum of 
the density profile of the nitrogen atom appears at a z distance 
less than 3A (Fig. |4](D)), indicating that molecules with the 
nitrogen atom pointing to the silica surface likely dominate at 



5 




■ N 

■ CT 
H 



10 15 20 25 30 35 40 45 

z(A) 
(A) 



5 10 15 20 25 30 35 40 45 

z(A) 




FIG. 4: Density profiles of the interfacial system defined from positions of all atoms treated as identical, positions of the center of mass, and the 
minimum position of acetonitrile molecules (A) and positions of the nitrogen atoms, alkyl carbon atoms, and hydrogen atoms of acetonitrile 
(B). The existence of three qualitatively different regions is evident: a LS region near the silica surface strongly perturbed by the substrate, a 
bulk liquid-like region from about 22 to 33A and a LV interface at larger z. C and D give an expanded view of A and B respectively from 
the surface to 6A . The density is normalized by the bulk value p B from a separate bulk simulation. The minimum position of an acetonitrile 
molecules refers to the minimum z position of all six atoms of the molecule. T=298K and CT stands for the alkyl carbon atom. 



close distances to the surface. 

We further divide the first layer into two sublayers accord- 
ing to the shoulder peak position of the density profile of iden- 
tical atoms z = 1.2A (Fig. 2](C)). A top view of the molecules 
within the sublayer < z < 1.2A was shown in Fig.|2](C) 
and (D). Clearly molecules in the first sublayer prefer to orient 
with their nitrogen atoms binding to hydrogen atoms (Fig. [2] 
(C)). About 25% of the sublayer molecules dock into the equi- 
lateral triangle formed by three closest oxygen atoms that have 
their attached hydrogen atoms pointing outside the triangle 
(Fig.[2](D)). 

Once we have defined the sublayer, the first layer, the bulk 
region, and the liquid-vapor region according to singlet den- 
sity profiles, it seems natural to examine the orientational dis- 
tributions in these regions as shown in Fig.[5](A). A consistent 
picture of the LS interface arises from all profiles. The orien- 
tation of an acetonitrile molecule in the first layer of the LS 
region tends to lie in two branches: 9 < 45° where the nitro- 
gen end points to the surface and around 9 ~ 130° where the 
nitrogen end points away from the surface. Most molecules 
within the sublayer z < 1.2 A lie in the first branch 9 < 45° 
while more molecules in the sublayer 1.2 A< z < 3 A lie 
in the second branch so the structure is reminiscent of a lipid 



bilayer. As the distance between an acetonitrile molecule and 
the silica surface increases, the orientational distribution grad- 
ually approaches unity, where the dipole of a molecule orients 
randomly as in the bulk. 

However, results in the LV interfacial region appear to de- 
pend on the choice of singlet density profile. When the loca- 
tion of a molecule is defined by its minimum position as in 
Fig. [5] (A), the orientational profile has a small maximum at 
9 = 90°, parallel to the interface. But Fig. [3(B) shows that 
the orientation profile of molecules defined using the positions 
of N, CT and cm seem to give different interpretations of the 
preferred molecular orientation at the LV interface. When the 
position of the center of mass is used, more molecules seem 
to orient parallel to the interface, consistent with results from 
the use of the minimum position of molecules in Fig. [5] (A). 
However, as noted in 11511 . more CH3 groups of the acetoni- 
trile molecules appear to point into the liquid (gas) phase if 
we define the location of an acetonitrile molecule using the 
position of its nitrogen (carbon) atoms. 

We show here that this apparent discrepancy arises from the 
different contributions of "over-counted" dipoles, arising from 
the different definitions of molecular position. These consid- 
erations become important only when the size of the molecule 
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FIG. 5: Orientational profiles of the acetonitrile dipoles relative to 
the surface normal (positive z direction) when z is the minimum po- 
sition of the acetonitrile molecules (A) and the orientational profiles 
of the LV interface when z is the position of nitrogen atom, methyl 
carbon atom and center of mass (B). The Jacobian factor sin 9 has 
been scaled out in all plots of angle distributions, so the bulk sys- 
tem has a uniform distribution. 9 is the angle between the dipole 
orientation (CT <— N) and the positive z direction. In our geometry 
9 — implies that the CN group points to the liquid phase at the LV 
interface and to the silica surface at the LS interface. 



is a finite fraction of the interface width, and when they are 
taken into account results from the different definitions can be 
related to each other. 



A molecular +/— dipole is defined as "over-counted" when 
its positive side (CT atom in the case of acetonitrile) falls into 
the region of interest (35 A< z < 42 A for the LV interface), 
while its negative side (the N atom in this case) is still located 
outside of the region, with a similar definition for an over- 
counted — /+ dipole with the negative side inside and positive 
side outside. The difference shown in Fig. (B) using def- 
initions from the CT atom and N atom can be expressed as 

5ovcr, + /-(#) - .9ovcr, -/+(#)■ 

The over-counted contribution of +/— dipoles can also be 
related to the density and angular distribution of the dipoles 



(0) 



: 1 +L/2cos9 



dz Po(z)f (z,9) 



zi—L/2 cos 9 



dz p+(z) 



for 0^9 < tt/2 



z 2 -L/2cos6 



(3) 



dzp {z)f (z,9) 



22+-L/2C 



dz p+{z) 



for tt/2 6 ^ 7T 



where L is the length of the dipole (distance between CT and 
N). Here po(z) is the density profile defined using the center 
of mass position and p+(z) is the profile when the location 
of the molecule is defined according to its atom with positive 
charge (these functions closely resemble each other and either 
can be used in the calculation). fo(z, 9)d9 is the probability 
for a given dipole with its center located at position z having 
orientational angle in [9, 9 + d9}. When +/— and Z\lz% are 
exchanged simultaneously, Eq. ([3]) then gives the contribution 
of over-counted — /+ dipoles. 

In order to check the validity of Eq. (0), we first calcu- 
late the average orientational distribution around z = z\ and 
z = Z2 as shown in Fig. [6] (A). Figure HtB) shows the com- 
parison between the calculation from Eq. (0) (red dashed line) 
when the average orientational distribution is used as the in- 
put and when the explicit histogram distribution of the over- 
counted molecules is used (solid black line). The overlap be- 
tween the two lines confirms the validity of Eq. ([3]). In fact, 
Eq. (0) can be accurately approximated in this case by simply 
taking fo(z, 9) = 1. The computed over-counted contribu- 
tions are barely changed. As shown in Fig.|6jB), Eq. (0 gives 
an excellent description of the simulation results. 

Clearly, because the number of over-counted molecules at 
the high density side z = Z\ of the LV interface cannot be 
compensated by those at the low density side z = z-2 and the 
magnitude of <7 ovcr + /_ {&) is non-negligible, the orientational 
profiles computed from different definitions of the position 
of a molecule in the interface give apparently contradictory 
predictions for the prefered orientation of the molecules. 

In the case of a bulk system, the over-counting from the 
two boundaries exactly cancels each other and g ovcl - + /-{&) = 
5ovcr, -/+($)■ In the case of the LS interface, because values 
of po(z) in the numerator of Eq. (0 is much smaller than the 
average density of the interfacial region (e.g. po(z = 4.5A) ~ 
0), the magnitude of g ovcr .+/_ is very small and makes almost 
no contribution to the overall orientational profile in the inter- 
facial region (e.g. the first layer orientational profile). There- 
fore, there is little discrepancy for cases of bulk and solid- 
liquid interfacial systems. 

In case of the LV interface, when the length of the dipole 
(e.g. L = 2.6A for an acetonitrile molecule) is comparable 
to the interfacial width (e.g. Z2 — z\ ~ 7A), <? vcr.-/+ is 
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FIG. 6: Orientational profile of molecules at the boundaries around 
Z\ — 35 A and Z2 = 42 A (A) and contribution from the over-counted 
+/— molecules to the orientational distribution (B). The black solid 
line (Simul.) in (B) is computed from the histogram of the over- 
counted molecules. The red dashed line (Theory) in (B) is computed 
through Eq. I0 using p+(z) ~ Po(z), shown in the inset of (B) (see 
also Fig.[4}, and f (z, 6) shown in (B) with z\ = 35A, z 2 = 42A 
and L = 2.6A. The green solid line (Random) in (B) is computed in 
the same way as the red dashed line except that fo (z,8) = 1. 



asymmetric and the over-counted molecules from the dense 
boundaries can make a significant contributions to the overall 
distribution and thus generate an apparent discrepancy. These 
considerations should be taken into account in interpretations 
of simulation or experimental data. 



B. Correlated structural properties 

Although singlet density and orientational profile have been 
studied in previous work flll[l3ll . correlated structures are sel- 
dom investigated. In order to shed light on how the antiparal- 
lel correlation l25ll between two acetonitrile molecule changes 
due to the effect of inhomogeneity, we define a quantity called 
the radial angular distribution function: 

1 Nc Nc i 

9cm{r, 0) = —J2J2' (S(r - rij )S(0 - By)) , (4) 



i=l 3 = 1 



where N c is the number of particles contained in the three dif- 
ferent regions defined above: the bulk region, LV region, and 



the first layer of the LS region. The g cm {r,6) is essentially 
the normalized histogram over independent pairs of molecu- 
lar dipoles in the three regions of interest; the prime on the 
sum indicates that terms with i = j are excluded. Here r»j 
is the distance between the center of mass of the i-th and j-th 
molecule and 9ij is the angle between the dipole orientations 
of the two molecules. The radial angular correlation function 
is normalized to 1 at large distance r and the integral over 9 
gives the usual center of mass radial distribution function. 

Fig.|7]shows that acetonitrile at the LV interface has an even 
stronger antiparallel correlation than in the bulk. In contrast, 
molecules in the first layer of the LS interfacial region com- 
pletely lose the preferred correlation at 9 = 180°. The sil- 
ica surface not only aligns the singlet molecular orientation 
such that it prefers to point toward the normal of the surface 
(Fig- EJA)), but it also significantly destroys the strong anti- 
parallel pair correlations seen for bulk liquid acetonitrile. At 
the free liquid-vapor interface, the singlet orientation differs 
only slightly from that of the random distribution in the bulk. 
The dominant intermolecular interactions cause LV interfa- 
cial molecules to preferentially orient in an antiparallel man- 
ner. But local interactions from the silica surface are strong 
enough to significantly modify the pair structure of molecules 
in the first layer. Since solid acetonitrile has even stronger 
anti-parallel correlations that would be similarly disrupted at 
the silica surface this, along with the general modification of 
bulk structure at a surface, could provide a possible explana- 
tion for the decrease of the melting point when acetonitrile is 
confined in silica pores. 



C. Singlet and collective rotational dynamics 

Both singlet and collective rotational dynamics are of inter- 
est as they can be measured from spectroscopic experiments 
10, [nj. The time correlation function (TCF) of the rotation of 
the symmetrical axis (CT-N) is defined as 



1 Nc 



(5) 



where v ; is the rotational axis of the i-th molecule and N c 
is the total number of molecules counted in particular region 
of interest. The collective reorientational behavior of the bulk 
acetonitrile system has been studied by Ladanyi and cowork- 
ers lO but few workers have addressed the collective dynam- 
ics of interfacial or confined systems. 

The collective polarizability II of a group of molecules is 
the sum of polarizabilities of each molecule. The standard 
point dipole/induced point dipole (DID) model expresses the 
individual polarizability 7f(j) for molecule j as a sum of a 
molecular component (gas-phase polarizability a(j)) and an 
interaction-induced component 



7r(j) = <*{j) + a(j) 



N 



w(k), 



(6) 
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FIG. 7: Radial angular distribution g C m{r, 9) at T = 298K for the first layer LS region (A), the bulk liquid region (B), the interfacial LV region 
(C), and the results after angular integration over 8 (D). The center of mass position is used for the definition of three different regions. 6 is the 
angle between the dipole axis of two molecules and r is the distance between the center of mass positions of the two molecules. 



where the dipole-dipole tensor Tjk between molecule j and fc 
can be written as 

7> = ^( 3 -- T )U 3fc C7) 

and f = r/r is the unit vector along the direction r. Eq. dp 
can be solved self-consistently or by matrix inversion 113411 . 
An efficient first-order approach (CC1) corresponds to replac- 
ing 7f(fc) on the right hand side of Eq, © by the gas-phase 
polarizability a(k) lO . 

The total collective polarizability II for a region of 7Y C 
molecules is 

n = ^7f(i) = n I+n 2 , (8) 

j 



where IIo and II2 are the isotropic and anisotropic polarizabil- 
ity of the region of N c molecules respectively. Different from 
the bulk system, the average of the xz or yz component of the 
anisotropic polarizability in the nonuniform system we study 
here is not zero. In order to focus on just the relaxation part 
of the anisotropic collective polarizability, we define the TCF 
C 2 (t) as 

c 2 (t) = (Tr[(n 2 (o) - <n 2 » • (n 2 (t) - <n 2 »]) (9) 

Analysis of C 2 (t)js crucial because it is directly relevant to 
the OKE signal 05143711 . The notation used in this subsec- 
tion is the same as the work by Hu et. al. l37ll . The TCFs 
of anisotropic collective polarizability for different regions of 
our system are primarily studied here using both the all-order 
and first order DID model. The gas phase polarizabilities of 
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FIG. 8: Singlet (A) and collective (B) reorientational correlations for 
acetonitrile in different regions. See Eqs. $5$ and ([9). The red dashed 
line in (A) is computed from a separate simulation of the bulk system. 
Dashed lines in (B) are the corresponding CC1 model calculations. 
Rotational time scales obtained by linear fit of the logarithm plot after 
t > 2ps are shown in Table U The definition of different regions is 
the same as in Fig. [5] 



TABLE I: Estimate of singlet and collective rotational time scales for 
acetonitrile in different regions." 

axis/model LS 1st layer Bulk region LV region 
CT-N 24(4) ps 3.4(2) ps 1.9(1) ps 
Collective 8.1(2) ps 1.4(1) ps 1.1(1) ps 
CC1 9.6(2) ps 1.5(1) ps 1.2(1) ps 

"Numbers in parentheses are uncertainties in the final digit. Random errors 
from the least-square fits of the data from 2ps to 4ps (Fig.[8]A) or 2ps to 3.5ps 
(Fig. [8]B) are much smaller than systematic errors that likely arising from 
this limited fitting range. The latter are crudely estimated here based on small 
variations of input data. 



the experiment work of Kittaka and coworkers ^ . Compar- 
ison between the all-order calculations and the first order ap- 
proximation in Fig. |8](B) shows that the effect of the higher 
orders is almost negligible in the bulk and in the LV region, 
consistent with the work of Elola and Ladanyi fl33f1 . 

Higher order effects are a little more important for the col- 
lective reorientation of molecules in the first layer close to the 
silica surface. The values of the rotational time scales deter- 
mined by fitting the long time behavior (2-4 ps for singlet ro- 
tations and 2-3.5 ps for collective rotations) to an exponential 
form are shown in Table J] The value of the time scale for the 
self -rotational correlation in the bulk and LV region are t-q = 
3.4ps and tlv = 1.9ps respectively. These values are different 
from those reported by Paul and Chandra lfl3tl . 10.75ps and 
5.75ps respectively, because they determined the orientational 
relaxation time scale from the time integral of the correlation 
function, and finite size and equilibration effects would be ex- 
pected to affect both these estimates in different ways. Both 
work agrees however that the reorientational motion has been 
significantly enhanced at the LV interface relative to the bulk 
phase. 

Both singlet and collective reorientation of the first layer 
of acetonitrile close to the silica surface are highly hindered. 
However, the collective reorientation in the first layer is much 
faster than the self-rotation of a molecule in the same re- 
gion. Fourkas and coworkers measured the collective reori- 
entation time scales in the bulk and in the pores as 1.66ps 
and 26ps respectively at T=290K and 1.42ps and 19ps respec- 
tively at T=309K 0]. We estimated the experimental values 
for T=298K, 1 .54ps and 22 ps by simply taking the average of 
the values at the higher and lower temperatures. 

Thus, our calculated bulk value 1 .39ps is about 10% smaller 
than the experimental value 1.54ps. Because experiments on 
acetonitrile in pores measure components arising from both 
bulk-like and surface molecules, and can be affected by other 
aspects of the confinement, our computed time scale for the 
collective reorientation of the first layer LS acetonitrile, 8. lps 
seems in reasonable qualitative agreement with the reported 
experimental value of about 22ps. However, a quantitative un- 
derstanding of the different time scales in the OKE experiment 
and test of the triexponential model used |01 would require a 
similar theoretical investigation of the OKE signals in silica 
pores, which is beyond the scope of this paper. 



D. Electrostatic potential and the long range contribution 



acetonitrile used here are an = 5.80 A 3 and a ± = 3.65 A 3 
rf36ll . The polarization effect from the atoms of the silica sur- 
face is not considered here. 

Fig- [U shows the singlet and collective reorientation corre- 
lations in a logarithmic plot. A separate calcuation of C(t) 
from the simulation of bulk acetonitrile (red dashed line in 
Fig.[8jA)) overlaps the corresponding C(t) for the bulk region 
in the inhomogeneous system. The rotational time scale in the 
region with a distance of larger than 25A from the surface is 
almost the same as in the bulk. This finding is consistent with 



In order to shed light on the importance of the correct treat- 
ment of the long-ranged electrostatics in this interfacial sys- 
tem, we analyze the total equilibrium charge density p* ot (r), 
comprised of the fixed charges in the silica substrate and the 
induced equilibrium mobile charge density p1 cc {r) in the ace- 
tonitrile fluid. The associated electrostatic potentials are given 
by Poisson's equation (in Gaussian units) as 

V tot (r) = J dr'pUr') " ■ ( 10 > 
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FIG. 9: Separation of the electrostatic interactions. (A) The potential from a point charge (cr = 0) and from a Gaussian charge distribution 
with width a = 4. (B) The bare charge density of the mobile acetonitrile molecules and the corresponding smooth charge density. (C) 
The comparison between the smooth charge densities of the total system and the mobile acetonitrile molecules. (D) The total force in z 
direction and its long-ranged part. The insets show the acetonitrile charge densities (B) and the total force (D) in the region of the LV interface. 
Smoothing eliminates the simulation noise still visible in the bare results. The blue lines in (B) and (C) are the same but on different scales. 
See Eqs. (B), (fTUl and (19). 



and 

V acc (r) = J dr' pl cc (r') ■ — L^. (11) 

In other work ll3ll l38tl . we have argued for the utility 
of splitting the basic Coulomb 1/r interaction appearing in 
Eqs. (fTOb and (fTTT i into short- and long-ranged parts based on 
a "smoothing length" a of order a characteristic nearest neigh- 
bor spacing: 

1 . , . , erfc(r/cr) erf(r/er) 

-=v (r)+v 1 (r) = K -!—L + v / ' ■ (12) 

r r r 

Here vi(r) is the functional form associated with the electro- 
static potential due to a unit Gaussian charge distribution of 
width cr, defined as 

^W = ^5^3 ejC p(-^)' (13) 



implying a v\(r) given by the convolution 

t*(r) = / dr' PG (r') • * = 5*0M . (14 ) 
J |r-r'| r 

As shown in Fig.|9jA), v\ (r) is slowly-varying in r-space over 
the length scale a. 

We showed in Ref. l39ll that for bulk acetonitrile the forces 
from the long-ranged parts of the Coulomb interactions es- 
sentially cancel when a is properly chosen to be 4A, and a 
very accurate description of local pair correlations arises from 
considering the simpler "strong-coupling" system with only 
truncated Coulomb interactions vo(r). Moreover, as argued in 
Ref. lHHHol], in nonuniform systems, where this force cancel- 
lation does not occur, the charge densities that best reflect the 
long-ranged electrostatics are not the bare densities Ptot( r ) 
and pl cc (r) appearing in Eqs. ( [Tol l and (fTTT i but rather the 
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Gaussian-smoothed charge densities 

p&(r) = J dr'pl t (r')p G (\r-r'\), (15) 

and 

^ e (r) = J dr'pl ce (v')p G (\T-r'\), (16) 

generated from Eqs. ( [Tol l. (fTTT i and ( TT4T > when only the long- 
ranged Coulomb component «i(|r — r'|) is used: 

| drV? ot (r')-«i(|r-r'|) 

= /^'^(r')-^- (17) 

and 



the bare charge density is much smaller, and the long-ranged 
((j = 4 A) charge density can contribute as much as 50% 
(see the inset of Fig. [9] (B)). Analysis of the force in the z 
direction (Fig.|9](D)) supports this same conclusion. There- 
fore, the long-ranged part of the electrostatics plays a rela- 
tively more important role in determining the structure and 
dynamics of acetonitrile in the LV interface than near the LS 
interface, where most aspects of the local structure is domi- 
nated by short-ranged forces. 

Fig. [9] (C) shows a comparison of the total smooth charge 
density (red) and the smooth charge density of the mobile ace- 
tonitrile molecule only (blue). The overlap of the two densi- 
ties at distances larger than 10 A indicates complete screening 
of the induced surface dipole beyond about 10 A, consistent 
with our finding in previous subsections that the dynamics and 
structure of acetonitrile in the region between 25 A and 30 A 
have reached the bulk limit. 



Vf ce (r,a) 



dr'^ cc (r')-«i(|r-r'|) 
dv' pZ^') ■ T-^—rr 



(18) 



for the total system and the mobile acetonitrile molecules re- 
spectively. 

Thus we are able to rationalize the contribution of the 
long-ranged part of the electrostatics through analysis of the 
smooth charge density p q ° ot (r) and the corresponding electro- 
static potential V* ot (r, a). One virtue of this perspective is 
that because of the integration over the smoothing length a, 
Ptot( r ) I s ver Y slowly varying along the interface in the x 
and y directions, unlike the bare charge density. As a result, 
V* ot (r, a) = V* ot (z, a) to a very good approximation 11381 14011 
and we can focus only on forces in the z-direction given by: 



F z (z,a) 



dV\ ot {z,(j) 
dz 



(19) 



A plot of p1° t (z), plcc( z ) an d Fz(z, a) is shown in Fig. [9] 
The magnitude of the acetonitrile charge density close to silica 
surface (0 < z < 10 A) is on the order of 0.01 e A~ 3 and it 
dramatically reduces to the order of 0.0001 e A~ 3 at the LV 
interface. The smoothed charge density in the LS region for 
the relevant a = 4.0 A is much smaller, with a magnitude of 
order of only 0.0002 e A~ 3 . 

This implies that the long-ranged (cr = 4 A) part of the 
charge density only contributes 4% to the total charge den- 
sity in the region of LS interface. Nevertheless, the long- 
ranged forces play an important role in properly describing 
the surface-induced dipole, as discussed in detail in Ref. [38| 
and illustrated below. However, in the LV interfacial region, 



IV. CONCLUSION 

We have investigated the structure and dynamics of the 
acetonitrile system on silica surface and at its LV interface. 
The change of melting point, the slowing down of dynamics 
for molecules close to silica surface and the tilting angle of 
the molecule at the LV interface were all explained in detail 
through our studies of structure and dynamics. The antiparal- 
lel correlations in the LS (LV) region was weaker (stronger) 
than in the bulk. Both singlet and collective reorientational 
time scales in the first layer close to the silica surface are ap- 
proximately 6 — 7 times slower than the corresponding time 
scales in the bulk. The collective reorientational time scale 
is about 3 times faster than that of the singlet reorientation. 
Moreover, we provide a general explanation of the ambigu- 
ity that arises in determining the orientational profile using 
interfaces defined with different choices of the atoms in a 
molecular system. Our current study of the acetonitrile sys- 
tem provides a general framework for future MD simulations 
of LS and LV interfacial systems. We plan to investigate other 
typical systems such as quadrupolar (carbon dioxide), dipolar 
(propionitrile) and associating fluids (hydrogen fluoride and 
water) on the silica surface and at their liquid- vapor interfaces. 
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